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ABSTRACT 

Two-dimensional simulations of mixing and fall back in non-rotating massive stars have been carried 
out using realistic initial models for the presupernova star and assuming standard spherically sym- 
metric explosions of 1.2 x lO^^erg. Stars of 15 and 25 M© with both primordial and solar composition 
were modeled. The zero metallicity supernova progenitors were compact blue stars and the amount of 
Ray leigh- Taylor induced mixing in them was greatly reduced compared with what was seen in the red 
supergiants with solar metallicity. The compact zero-metal stars also experienced more fallback than 
their solar metallicity counterparts. As a result, the ejected nucleosynthesis from the two populations 
was very different. For the simple explosion model assumed, low metallicity stars ejected ejected 
too little iron and intermediate mass elements even to explain the abundance patterns in the most 
iron-poor stars found to date, suggesting that some important ingredient is missing. Rotation is likely 
to alter these conclusions by producing a greater fraction of red supergiants among Population III 
stars. The velocities of the heavy elements in all models considered - both red and blue supergiants - 
were less than observed in SN 1987A, suggesting that at least occasionally, asymmetric aspects of the 
explosion mechanism and fallback play a major role in mixing. 
Subject headings: Supernovae, nucleosynthesis, first stars 



1. INTRODUCTION 

The nucleosynthetic yields of supernovae are important 
components of galactic chemical evolution and are essen- 
tial to understanding the abundances observed in metal- 
poor stars. The yield of a core-collapse supernova is 
determined by its presupernova evolution, the geometry 
and energy of its explosion, and by the mixing and fall- 
back that occur as the supernova shock traverses the star 
(e.g. , I Chevalier 1 2005 1 ) . Mixing and fallback also affect the 
supernova's light curve and spectrum and the appearance 
of its re mnant. The dense knots and fila ments visible in 
Hubble ("Blair et al. 2000) and Chandra (Hughes et al. 
jpOOO) images of supernova remnants provide clear evi- 
dence that mixing of some sort is a common occurrence. 
Supernova 1987A, the closest and most thoroughly ob- 
served modern supernova, provided important observa- 
tional constraints on models. Its smooth bolometric light 
curve required extensive mix ing of the heliu m core with 
the hydrogen envelope (e.g., Woosley |1988). The early 
appearance of X-rays and 7-rays from radioactive decay 
of ^^Ni, as well as spectroscopic evidence some fraction of 
the iron peak was mixed out to 4000 km s~^, provided 
further evidence for extensiv e mixing in the interior of 
the supernova at early times ( Arnett et al.||1989 Witte- 
[born et al. 1989; Utrobin 2004^ : 

Differences m the presupernova structure of zero and 
solar metallicity stars alter the way mixing and fallback 
processes operate. Below about 40 M©, non-rotating, 
zero-metallicity stars are expected t o be more com- 
pact t han solar-meta llicity progenitors ( Heger & Woosley 



nuclei for the CNO cycle, leading to inefficient hydro- 
gen burning and a very dense hy drogen shel l of low en- 



20081 [Hirschi et aT] [2008 ) . Metal-free gas has a lower 
opacity than solar-metallicity gas and lacks initial seed 
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tropy in the presupernova stars. [Chevalier (1989) pre- 
dicted that higher amounts of fallback are expec ted for 
more compact progenitors. In a recent paper, | Zhang[ 
[et al.| (2008) used their one-dimensional Eulerian code 
PANGU to determine the remnant masses left behind by 
the sup ernova models calculated in su rveys by j Woosley &[ 
Hegerf([2QQ7) and Heger & Woosley[ ([2008[). They found 



that zero-metallicity supernovae experienced more fall- 
back and left behind larger compact remnants than their 
solar metallicity counterparts. For example, the baryonic 
remnant masses left behind by 25 Mq stars of zero and 
solar metallicity were 4.16 and 2.09 M©, respectively. 

Different presupernova structures, arising from differ- 
ences in stellar mass and metallicity, determine where 
Ray leigh- Taylor instabilities occur and the extent to 
which they grow. An initially static, incompressible fluid 
is unstable if the pressure gradient points opposite to 
the density gradient, i.e., when (dP/dr)(d p/dr) < (e.g. 
[Chevalier'1976! [Benz fc Thi elemann'1990). The location 
oTthese density mversions varies with time as the forward 
and reverse shocks propagate through the star. Partic- 
ularly important are regions where the forward shock 
encounters an increasing value for pr^, where p is the 
density and r, the radius (Herant & Woosley" 1994"). The 



time scale also depends upon the initial stellar structure. 
In particular, a more compact star will experience faster 
shock propagation, leaving less time for instabilities to 
grow. 

Because SN 1987A was so well observed, most previous 



novae ([Arnett et al. 1989| 


Fryxeh et al.|[1991| Mueller 


et al.||1991| |Hachisu et al.| 
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, 1992| [Herant & Benz 
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Kifonidis et al.[ 


2006 


) have been in the con- 



progenitors ( Herant fc Woosley|l994 ). Except for Herant 
and Woosley, these studies all used red and blue super- 
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giants of 15 to 20 Mq as progenitors. The methodol- 
ogy of ah these studies was similar. A one-dimensional 
progenitor model was exploded, somewhat artificially, by 
means of a piston or a thermal bomb, and the subsequent 
evolution followed with a two-dimensional code. 



More rece ntly, [Kifonidis et al., (2003) and Kifonidis 
et al. (2006) have used a different approach. These au- 
thors followed a blue supergiant model from the first sec- 
onds of the explosion out to about 5 days after core col- 
lapse, using first one code with neutrino physics for the 
early times, and another code with mesh refinement for 
later times. They saw mixing at the Si-0 interface, a 
location at which no previous studies had found mixing, 
and were able to reproduce the high ^^Ni velocities ob- 
served in SN 1987A, something previous studies had not 
done. 

While attempts to reproduce observations of 1987A 
have been numerous, no multidimensional studies of mix- 
ing and fallback in very low metallicity supernovae have 
been done. The nucleosynthetic yields of metal-free (Pop 
III) and extremely metal-poor (EMP) stars might still be 
visible in the abundance patterns observed in some halo 
stars in our own galaxy. Of p articular interest are the 
"ultra-iron- poor' 
2006^. 



et al. 



(HMP) stars ( |Frebel et al 
These stars with [Fe/HJ< 
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lave abun- 



dance patterns that differ considerably from those ob- 
served in stars wi th near solar meta llicity or even other 
metal-poor stars ( Cayrel et al.||20Q4| ). It is possible that 
these iron-p oor stars were enric hed by only one or a few 
supernovae (Frebel et al. 2005). In particular, the two 
most metal-poor stars known and several other UMP 
stars display marked enhance ment in C, N, and O rel- 
ative to Fe. Previous studies (|Iwanioto et al.|2005 Tom 



inaga et al.|2007[ |Heger fc WoQsIey|2008p have sought to 
explain these abundance patterns with one-dimensional 
models for supernovae that parametrize the amount of 
mixing and fallback to match what is observed. Simulat- 
ing mixing and fallback directly, rather than parametri- 
cally, requires a mult i- dimensional approach. 

In this paper, we use two-dimensional axisymmetric 
simulations to explore directly the amount of Rayleigh- 
Tay lor- induced mixing that occurs in non-rotating zero- 
and solar- metallicity stars. Our methodology is similar 
to the earlier studies (before Kifonidis) of SN 1987A. In 
§[2j we discuss our initial models, our modifications to the 
FLASH code, and our simulation setup. In §[3) results are 
given that show the degree of mixing, the final velocity 
distribution of isotopes, and the ejected yields. These 
yields are compared with abundances observed in HMP 
stars in § |3.4[ Finally, we provide a short summary of 
results and their interpretation in § |4] 

2. MODELS AND METHODS 

The present work follows the method used in many pre- 
vious studies of mixing in supernovae. A one-dimensional 
code was used to evolve and explode the pre-supernova 
model and to follow the first stages of the expansion 
to the time when the reverse shock was just beginning 
to form. No significan growth of instabilities is ex- 
pected before the formation of the reverse shock. The 
one-dimensional model was then mapped onto a two- 
dimensional grid and the ensuing instabilities followed. 

While Kifonidis et al.| ([2006 ) reproduced the observa- 
tions of SN 1987 A somewhat better than previous at- 



tempts, possibly by following the early stages of the ex- 
plosion, this paper does not do that. The physics of the 
initial explosion remains uncertain and our goal is to iso- 
late the differences in post-explosive mixing that arise as 
a direct consequence of the differences in initial structure 
of the pre-supernova models. Exploding the star with a 
piston in the same location and with the same energy in 
all models allows us to accomplish this. 

The choice of piston mass location is constrained by 
observational parameters. The piston cannot be located 
within the iron core or the resulting explosion will pro- 
duce far too much of 54,58 Fe and other neutron-rich 
species to be in agreement with observations of these 
isotopes. On the other hand, the remnant mass will be 
too large to agree with observations if the explosion site 
is located outside the base of the oxygen shell. There 
are reasons to believe the location site is located at the 
base of the oxygen shell-the large density decrease asso- 
ciated with this location is dynamically important, and 
successful explosion calculations often find the mass cut 
there. 

The papers from which our models are taken also re- 
ported on models in which the stars were exploded with 
a piston at the edge of the ni core. These explosions ex- 
perience slightly less fallback, and produce more nickel, 
by a factor of 2 or so, than the models presented in this 
paper. 

2.1. Progenitor Models 



Initial models were taken from the surveys of Heger 
Woosleyl ( [2008) and Woosley fc Heger ( 2007). Both~ 



these papers used th e KEPLER code ( [Weaver et al.||1978 
Woosley et al. 2002) to evolve stars through all stable 



stages of nuclear burning until their iron cores became 
unstable to collapse. At this point, pistons located at or 
near the bas e of the oxygen shel l were used to explode the 
stars. Heger fc Woosley ( 2008| simulated the evolution 
and explosion of 10 to 100 Mq stars with zero initial 
metallicity. Explosion energies range d from 0.3 to 10 
where 1 Bethe = 1 B = 10^^ ergs. Woosley fc Heger 



(2007) examined solar- metallicity stars from 12 to 100 
M0 which were exploded by pistons similar to the other 
survey, but for a more limited set of masses and energies. 
Both surveys were limited to non-rotating progenitors. 
The zero-metallicity stars were assumed to have no mass 
loss, while the solar-metallicity models took mass loss 
into account. 

Here only two representative stars fr om each survey are 
studied: Models zl5D and z25D from [Heger fc Woosley 



(2008) and Models sl5A and s25A from'Woosley fc Heger 
(j2007). The letter "z" indicates zero initial metallicity. 
while "s" indicates solar metallicity. The numbers in the 
models correspond to the initial mass of the star in M© 
and the final letter is the explosion energy, 1.2 B in each 
case. The piston was located at the place in the star 
where the entropy was equal to 4.0/cB/baryon. This cor- 
responded to the base of the oxygen shell. Series sA and 
zD are thus directly comparable in all respects save their 
initial metallicity. 15 and 25 Mq represent the "canoni- 
cal" supernova cases with the most commonly employed 
explosion energy and piston location. 15 and 25 Mq 
stars are also in the same mass range as previous stud- 
ies of Rayleigh-Taylor mixing in supernovae, allowing for 
easy comparison. 
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The one-dimensional progenitor models used in this 
study lead to spherically symmetric explosions. Map- 
ping the models from one to two dimensions after the 
explosion has taken place has the effect of suppress- 
ing low-order departures from spherical symmetry. This 
work only addresses spherically symmetric explosions 
with asymmetries of significantly higher mode than / = 1 
or 2. 

2.2. The FLASH code 



The FLASH code (iFryxeh et aT]|2QQQ! ) was used to fol- 
low shock propagation and mixing m the models in two 
dimensio ns. This code has been extensively ver ified and 
tested ( Calder et aI]|2QQ2[ [Weirs et al^|2QQ5[ ). FLASH 
is an adaptive mesh refinement code based on an Eule- 
rian implementa tion of the Piecewise Parabolic Method 
(PPM) of Col ella fc Wood ward (1984). The code can be 
configured in a number of different ways. We used the 
HLLE Riemann solver to resolve shocks. KEPLER uses 19 
different isotopes to evolve stellar models through stable 
stages of nuclear burning. We use the "aproxl9" compo- 
sition module included with the FLASH2.5 distribution to 
map these isotopes directly to FLASH. The abundances of 
19 different isotopes, from to ^^Ni, are stored for each 
grid cell. Explosive nuclear burning could have been fol- 
lowed in FLASH using this network, but the burning was 
over by the time the models were mapped into FLASH. 
Mapping the star in at earlier times to follow the ex- 
plosive burning in two dimensions would have had little 
effect, as no departure from spherical symmetry is ex- 
pected until the formation of the reverse shock, which 
occurs much later in the calculation. FLASH was con- 
figured to use axisymmetric coordinates. The gravita- 
tional potential was computed using a multipole method 
to solve Poisson's equation. The mass distribution in 
this simulation had only slight deviations from spheri- 
cal symmetry, so the additional gravitational force from 
overdensities in the simulations was small. The gravi- 
tational potential from the point mass at the origin of 
the grid was added to the potential computed by the 
multipole solver at each time step. The hydrodynamic 
equations were solved using an explicit, dimensionally- 
split approach: the hydrodynamic equations were solved 
first along one coordinate grid direction, then the other 
at each time step. 

An equation of state was employed that assumed full 
ionization and included contributions from radiation and 
ideal gas pressure: 
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where P is the pressure, a is the radiation constant, 
is Boltzmann's constant, T is the temperature, p is the 
density, rup is the proton mass, p is the mean molecular 
weight, and E is the energy. Although the outer regions 
of the remnant may not always be fully ionized, the re- 
gions where Ray leigh- Taylor mixing takes place are. Ex- 
ploratory simulations performed in one dimension with 
FLASH using a Helmholtz equation of state were identical 



to one-dimensional simulations performed with the per- 
fect gas and radiation equation of state used in the sim- 
ulations presented in this paper. These one-dimensional 
simulations were only run to 2 x 10^ seconds. At later 
times, the density is low enough that the Helmholtz equa- 
tion of state no longer applies. 

The code was configured to use pressure, density, ^He, 
and ^^O as its refinement variables. An error estimate for 
a block was computed using the second derivative of the 
chosen refinement variables. If the estimated normalized 
error in one or more of these variables was greater than 
a given value, regions were refined until a normalized 
error is reached that is less than the acceptable value 
or the maximum level of refinement is reached. Regions 
of the simulations with steep gradients in one or more of 
these variables were likely to be refined. If the normalized 
error estimate was below a certain value, i.e., the absolute 
value of the second derivative of one or more variables 
was small, then the region was "de-refined" . 

2.3. Modifications to the FLASH Code 

The FLASH2.5 release was customized to include a mod- 
ule that inserts a roughly circular zero-gradient inner 
boundary around the origin. This prevented infalling 
matter near the origin in the simulations from backing 
up and affecting the outer regions. Matter was allowed to 
fall though the zero-gradient, quasi-circular boundary at 
the center of the model and accumulate on the point mass 
at the origin. If the radius for the inner boundary passed 
through one of the inner zones, the zones interior to that 
zone were set to be duplicates of the zone on the bound- 
ary. While the Cartesian nature of the axisymmetric co- 
ordinate system meant that the inner boundary was only 
as close to circular as one can reproduce with square com- 
ponents, it did not introduce a significant amount of error 
into the calculation. The inner boundaries were chosen 
to be within the sonic radius, ensuring that small numer- 
ical errors at the interior boundary did not accumulate 
and affect the flow of fluid upstream from the boundary. 
The sonic point moved outward, not inward, as the simu- 
lation time progressed, so the inner boundary was always 
within the sonic radius. The falling temperature caused 
the sound speed to decline with time, while the velocity 
near the inner boundary increased with time. The gravi- 
tational potential resulting from this point mass was up- 
dated and added to the gravitational forces computed by 
the multipole solver at each time step. 

A module was also included that locally deposited en- 
ergy from radioactive decay of ^^Ni to ^^Fe. The amount 
of energy 

d£;56Ni = A56NiX56Nie-^^^Ni*^ * q{^^m) (3) 

The decay rate of ^^Ni, Ase^i, is 1.315 x 10"^ s"^, and the 
amount of energy released per gram of decaying ^^Ni is 
g'(^^Ni), for which we took the value 2.96 x 10^^ erg g~^. 
X56Ni is the fraction of ^^Ni in the block. The amount of 
^^Co at a given time could be found as a function of the 
amount of initial ^^Ni by 
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co-t) (4) 



^56Co — Ase^i 

SO that the energy deposition rate from ^^Co as a function 
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of time was given by 



dE56Co 



■>Co - A56N: 



-X56Ni(e-^^^Ni^_e-^^6co^))A56Co*g( 

(5) 

We assumed a decay rate for ^^Co, Ase^o, of 1.042 x 10 
s~^, and an energy per gram of decaying ^^Co, q{^^Co)^ 
for which we took the value of 6.4 x 10^^ erg g~^. 

2.4. Calculations 

The supernova models were evolved with KEPLER to the 
point where all explosive nuclear burning had ceased and 
the reverse shock had just begun to form. This occurred 
at 10^ s and 10^ s for 15 and 25 stars of solar com- 
position, and at 25 s and 100 s for 15 and 25 Mq stars 
of primordial composition, respectively. At these times, 
the one dimensional models from KEPLER were mapped 
onto a two-dimensional axisymmetric grid and evolved 
forward in time with the FLASH code. A similar simula- 
tion of Model s25 A was also performed using a progenitor 
evolved to 10^ seconds with KEPLER before being mapped 
to FLASH. No significant difference in the later evolution 
of s25A models evolved to 10^ seconds and 10^ seconds 
was observed. Only one quadrant of the star was carried 
in the calculation, enforcing symmetry about the left y- 
and bottom x-axes, while allowing material to leave the 
grid through a zero-gradient boundary at the right y- 
and top x-axes. 

An enhanced flow was observed along the x- and y- 
axes. This flow was not large in comparison with the 
rest of the simulation, but it was present, as can be seen 
in Figures [7| - [6j This is a well-documented artifact of 
the dimensionally-split approach to solving the hydrody- 
namic equations. It did not substantially influence the 
evolution of the simulations presented here. 

Perturbations arising from a Cartesian grid are also 
inevitable. In order to quantify these grid effects, we 
performed simulations of all stars with a random pertur- 
bation in velocity with a maximum amplitude of 0.5% 
and 2.0%. We also performed simulations with no ad- 
ditional perturbation. We found that perturbations of 
2.0% in velocity had a clear effect on the initial scale of 
the Ray leigh- Taylor instability, increasing the amplitude 
and scale of the flrst instabilities to form. This effect 
is shown for si 5 A in Figure [ij The case of 0.5% ran- 
dom perturbations to the velocity results in a scale for 
the initial Rayleigh- Taylor instabilities that is roughly 
equivalent to the case where the only perturbations were 
those arising from the Cartesian grid itself. Other models 
had similar resolution, such that grid perturbations were 
roughly equivalent to velocity perturbations of 0.5% and 
velocity perturbation of 2.0% had a noticeable effect on 
the initial scale of the instability. Perturbations arising 
from the grid were no larger than 2%, well within the 
regime of perturbations expected to arise from convec- 
tion. 

Because the zD-series models were more compact than 
sA-series models, the reverse shock reached the centers 
of zD models faster than in the sA models. This had 
the effect of shutting off mixing in the zD-series before 
the Rayleigh- Taylor in stability had time t o become fully 
non-linear (see also Her an t fc Woosley[ 1994). The initial 
perturbations had a greater effect on the final state of the 
zD simulations, while the initial perturbation scale and 
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spectrum are washed out in the sA models as a result 
of their longer mixing times. A simulation with random 
C^locity perturbations of 5% was performed for Model 
z25D. 

Models were initially mapped onto the two- 
dimensional grid such that their inner iron cores 
were resolved with at least 4 blocks of 16 zones each. 
This was sufficient to ensure that the rest of the star was 
accurately resolved. As the simulation progressed and 
the model stars expanded, the maximum refinement level 
of the simulation was turned down, so that the model 
star was always resolved at about the same percentage of 
the radius of its inner core of "^He and heavier isotopes. 
This is the region where Rayleigh- Taylor mixing takes 
place. The star expands homologously, ensuring that all 
regions will be adequately resolved. 

The solar composition models were mapped onto a grid 
5 X 10^^ cm on a side. The zero-metallicity stars were 
mapped onto a grid 1.4 x 10^^ cm on a side. The portion 
of the grid outside the original KEPLER model was initial- 
ized with a density proportional to r~^. Simulations with 
an outer density proportional to r~^-^ was also performed 
for Models sl5A, s25A, and z25D. These density profiles 
span the realistic range of smooth density distributions 
outside real stars. No difference in the final profiles for 
density, temperature, pressure, or composition was found 
between models with different outer density profiles. The 
density profile of the surrounding material therefor has 
no effect on the amount of Rayleigh- Taylor mixing that 
goes on inside the star, provided very little mass as a 
proportion of the original mass of the star is swept up 
in the first days of the explosion. The amount of mass 
added to the grid from ambient density was < 2% for all 
models. 

Calculations were run at least until the Rayleigh- 
Taylor fingers had ceased to move with respect to the 
mass coordinate of the star. This happened 2 hours after 
core bounce for Model zl5D, 4 hours after core bounce 
for z25D, and 7 days after core bounce for the sl5A and 
s25A models. All models were followed to 10^ seconds, 
long after Rayleigh-Taylor mixing had frozen out in the 
zero metallicity stars, but long enough that infall though 
the inner boundary had reached an asymptotic stage and 
the final remnant mass from these two-dimensional sim- 
ulations could be determined. 

3. RESULTS 

Figure [2] shows the evolution of the stable and unstable 
regions of the models, and the position in mass coordi- 
nate of the forward and reverse shocks. A reverse shock 
forms when the outgoi ng shock en counters a regio n of in- 
creasing pr^ ( [Herant fc Woosley 1994[|Woosle y fc We aver] 
,1995, ) . When the shock encounters a density profile that 
falls off with a flatter slope than ex r~^, i.e. when it en- 
counters a region of increasing pr^ it decelerates. The 
deceleration of the forward shock front reverses the di- 
rection of the pressure gradient, which slows down the 
layers interior to the shock, as well. Shocked material 
piles up and forms a high density post-shock shell. The 
reverse shock forms at the inner boundary of the high- 
density shell of decelerated matter and propagates down 
into the star, toward its center, slowing down the deeper, 
inner layers of the star (Kifonidis et al. 12003). The decel- 



eration of the shock creates a steep pressure gradient in 
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Fig. 1. — Influence of perturbations on the initial scale of the instability. Panel (a): no additional perturbation. Panel (b): random 
perturbations of 0.5% applied to velocity. Panel (c): random perturbations of 2.0% applied to velocity. Panel (a) and (b) show instability 
growth on nearly the same scale, implying that grid perturbations were on the order of 0.5%. The scale of the instability in panel (c) is 
noticeably larger, showing that perturbations from the grid are no larger than 2.0%. The other models presented in the current paper, not 
shown, are nearly the same. The initial scale and spectrum of the instability is washed out after the instability becomes nonlinear in the 
solar stars. 
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Fig. 2. — Stability evolution for the flrst ~ 5% of simulation time. 
Left axes show the location of the base of the shell indicated. Note 
different time and mass axes scale between the models. The reverse 
shock takes far longer to propagate back for solar composition than 
for Pop III stars. A greater portion of mass in the solar stars is RT 
unstable, leading to more mixing in these stars. 

the opposite direction to the existing gravitational and 
density gradients. This pressure gradient can overwhelm 
the gravitational gradient, and in doing so triggers the 
formation of Rayleigh-Taylor instabilities in the mate- 
rial. The Rayleigh-Taylor instability develops until the 
reverse shock has passed by, at which point the material 
becomes stable again, and the instabilities cease to grow 
exponentially. Figure |2] covers the period of time from 
when the models were first mapped to FLASH to slightly 
beyond the time when the Rayleigh-Taylor instability fin- 
gers began to grow. 

Figure [2] shows that the Rayleigh-Taylor instability 
had far more time to develop in the solar composition 
models than primordial composition models. The zero- 
metallicity models, particularly Model zl5D, were far 
more compact than their solar metallicity counterparts. 
The reverse shock took only ^100 seconds to reach the 
center of Model zl5D, leaving very little time for the 
Rayleigh-Taylor instability to grow. This is reflected in 
the low degree of mixing seen in Figure [4j which shows 
the final distribution of isotopes as a function of mass 
for all stars. The original structure of these stars can be 
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Fig. 3. — Distribution of isotopes just after core bounce as a 
function of mass coordinate. The dotted vertical line shows the 
position of the mass cut-matter to the left of this line will fall 
back onto the star. Without mixing, no ^^Ni escapes from stars 
of primordial composition, while most of the ^^Ni core is expelled 
from the solar composition supernovae. 
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Fig. 4. — Distribution of isotopes after mixing has ceased. Note 
that we flnd almost no difference between this distribution and the 
initial conflguration for Model zl5D. Mixing is conflned to the O- 
He shells for Model z25D. Mixing has penetrated the Si/O layer 
for Model sl5A. 
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Fig. 5. — Log(density) snapshot of Model zl5D after mixing 
has ceased. The RT instabihty visible in the upper right hand 
corner is most likely a result of the artificial outer density profile. 
It represents the outer layer of the star. The mixed inner layer 
extends only about 1/10 of the radius of the star, covering a much 
smaller region than in the solar models. 

seen in Figure [Sj Two-dimensional snapshots of mixing 
are shown in Figure [5j which shows the density struc- 
ture of the entire star, and Figure |9) which shows the 
isotopic composition of the mixed region at the center of 
the model. In Model z25D, the reverse shock took longer, 
about 10^ seconds, to reach the center of the star, allow- 
ing the RT instability to grow for a longer period of time. 
Figured shows an unstable band between the ^He/^^C- 
^^O shell boundary, and that was indeed the place we 
saw mixing in these models. More mixing in z25D than 
zl5D can be seen in Figures |4j [6j and[lOj 

Solar composition models were about 50 times larger 
in radius than their primordial composition counterparts, 
and correspondingly less dense. The reverse shock took 
longer to form and 10^ and 2 x 10^ seconds to propagate 
back to the mass-coordinate origin for Models sl5A and 
s25A, respectively. This was about 100 times longer than 
for Model z25D, giving the Rayleigh-Taylor instability 
more time to develop. Additionally, a wider range of re- 
gions between and including the ^He and Si+S shells were 
unstable than in the zero metallicity counterparts. Fig- 
ure [2] implies that the solar metallicity stars were mixed 
to a greater degree than stars of zero metal initial com- 
position. Figures [7| and |8] show that about half the solar 
composition stars were mixed, compared to only about 
1/10 of the zero metallicity stars. 

3.1. Infall 

Infall is as important as mixing for determining the 
final yield of a supernova. For a freely expanding super- 
nova remnant, the rate at which mass accretes onto the 
black hole or neutron s tar at the center of the explosion 
IS given by M (X (|Chevalier||l989|. It takes super- 

nova explosions on the order of 10° seconds to reach this 
asymptotic, freely expanding stage, at which time it is 
possible to determine the final mass of the stellar rem- 




FlG. 6. — Log(density) snapshot of Model z25D after mixing has 
ceased. Like Model zl5D, a small portion of the star is mixed when 
compared to the corresponding model of solar composition. The 
RT instabilities visible in the upper corner formed at the boundary 
of the star with an artificial density background. 
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Fig. 7. — Log(density) snapshot of Model z25D after mixing has 
ceased. The Rayleigh-Taylor instability that began at the O/He 
interface has mixed the inner layers of the star down to the iron 
core. The fiow visible along the y-axis is a result of dimensional 
splitting in the hydrodynamic solver. 

nant by extrapolating from the asymptotic infall rate. 
Our final remnant masses are compared with those found 
with one-dim ensional calculations carried out by , Zhang] 



et al.| ( |2008| ) with the PANGU code in Figure [13 

PANGU IS a one-dimensional hydrodynamics code based 
on the second-order semi-discrete finite difference cen- 



tral scheme of Kurganov & Tadmor (2000). Time evo- 
lution was carried out by a third- order total variation 
diminishing Runge-Kutta method (Shu & Osher 1989). 
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Fig. 8. — Log(density) snapshot of Model s25A after mixing has 
essentially ceased. The Ray leigh- Taylor instability that began at 
the O/He interface has not penetrated through the thick ^^0[ shell, 
as it had with model sl5, where the shell was thinner (Figure [t]) 




Fig. 9. — Final distribution of isotopes in the inner core of 
Model zl5D. Note that the mixed region extends only to about 
1/10 of the radius of the star-a much smaller proportion than in 
the solar metallicity stars. Pictured is the run with a 2% percent 
perturbation in velocity. The Ray leigh- Taylor instabilities have 
had a shorter time to grow, and so have retained more of their 
original shape, than in the solar metallicity models. 



Zhan^ et al. (2008) simulated the explosi ons of stars in 
the |Woosleyl^ Heger,(2007^ and ,Heger fc Woosley| ( [2QQ8 ) 
surveys. Using this one- dimensional code, the authors 
followed the evolution of the supernova remnant out to 
10^ seconds, to the time at which the accretion rate onto 
the central remnant had reached an asymptotic depen- 
dence on time and the final remnant mass could be de- 
termined. 






Fig. 10. — Final distribution of isotopes in the inner core of 
model z25D with a 5% perturbation in velocity. Note the mixed 
region goes out to only about 1/10 of the radius of the star-a much 
smaller proportion than in the solar metallicity stars, although the 
mixed region is greater and the initial scale of the perturbation is 
less visible than in Model zl5D. The reverse shock took longer to 
pass though Model z25D then Model zl5D. 






(d) 



Fig. 11. — Final distribution of isotopes in the inner core of Model 
sl5A. The helium, oxygen, and silicon shells have been disrupted 
by the Rayleigh-Taylor instability, which has mixed ^^Ni out of 
the core and -"^H all the way in to the inner layers. Model sl5A 
is more completely mixed than any other model presented in this 
paper, because of both the long timescale for the reverse shock and 
a thinner oxygen layer. 

Our final remnant masses showed good agreement with 
the PANGU results for all Models except Model z25D, as 
shown in Figur€ 13 The shape of the infall curves for our 
stars, as shown in Figure [Tsj match the PANGU results 
well, though the remnant masses from FLASH are larger 
except in the case of s25A. Model s25A was mapped to 
FLASH after 1 x 10"^ seconds of evolution in KEPLER. Al- 
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Fig. 12. — Final distribution of isotopes in the inner portion of 
Model s25A. The ^He layer has been completely disrupted, though 
the RT instability has not penetrated past the ■'^^O layer. The ring 
of ^^Ni visible was formed by explosive nuclear burning. ^^Ni has 
not been mixed out of the center of the star, and the -^^Si layer is 
only marginally effected by the Ray leigh- Taylor instability. 

though this does not alter the final mixed state of the 
star, it does alter the amount of mass that accumulates 
at the inner boundary, since most of the infall occurs 
during the first 10^ seconds, when the star was still be- 
ing evolved forward with the Lagrangian code KEPLER. 
Both 15 solar mass models are in good agreement with 
the PANGU results. z25D shows the most deviation, most 
likely because this model experienced the most fall back 
of any of the models studied here. We were unable to 
replicate the fine resolution employed by a one dimen- 
sional code, and that is almost certainly the reason for 
the enhanced infall mass. The initial remnant mass in 2D 
was sufficiently larger than the ID remnant mass that the 
additional force it exerted on the surrounding material 
was large enough to cause additional infall. This led to 
a larger remnant, causing a mild runaway effect. While 
the remnant mass obtained with FLASH for Model z25D is 
likely inaccurate, Rayleigh-Taylor mixing in this star has 
stopped by 2 X 10^ seconds, at which point the remnant 
mass in the FLASH simulation was not significantly larger 
than that in the PANGU simulation. It is unlikely that 
increased infall has had a significant impact on the evo- 
lution of the Rayleigh-Taylor instability, so the mixing 
results remain sound. 

3.2. Perturbations 

Random perturbations of amplitude 0.5% and 2% of 
the original velocity profile were added to the initial 
mode l for the solar composition models, as described in 
Random perturbations of 2% were added to the 



2.2 



zero-metal stars, and an additional simulation with an 
initial perturbation of 5% was performed for Model z25D. 
Perturbations of 0.5% produce an original size scale for 
the Rayleigh-Taylor instabilities that is about the same 
size as those arising for the non-perturbed case, while 
the 2% perturbations result in a larger initial scale for 




Fig. 13. — Comparison of infall rates through the inner boundary 
for 1-d PANGU and 2-d FLASH Eulerian codes. For all stars except 
Model sl5A, FLASH overestimates the initial infall rate, which is 
reflected in a larger remnant mass. This is due to limited resolution 
at the origin in the FLASH code. The shape of the FLASH infall 
curves matches the shape of those curves in PANGU. The exception 
is Model z25D, which shows a dramatic departure from the shape of 
the PANGU curve at later times. This is probably because the extra 
remnant mass is large enough to effect the infall rate through the 
inner boundary. 




Fig. 14. — Mass fraction of chemical isotopes as a function of 
radial velocity for Model zl5D for different perturbations. Larger 
initial perturbations lead to wider distributions. These differences 
effect about 10% of a given isotope. 

the instability, implying that grid perturbations for these 
models are between 0.5% and 2%, as shown in Figure [l] 
The Rayleigh-Taylor instabilities can grow for at least as 
long as the reverse shock takes to reach the origin. The 
instabilities in the solar metallicity stars can grow for 
many e— folding times, a long enough time to wash out 
the initial scale and spectrum of the instabilities. The 
final states of the perturbed models appear essentially 
identical to the unperturbed models. The distribution of 
the isotopes in both velocity (Fig ures 16 and 17) and and 
mass space (Figures 20 and 21) has no systematic cor- 
relation with the magnitude or the initial perturbation 
for solar metallicity stars. The scale of the perturbations 
one would expect in a real star is set by the magnitude of 
the convective velocities, which are on the order of 0.5% 
of the total velocity. 

Perturbations for the zero-metallicity case have a 
greater effect on the final amount of mixing in these stars. 



9 



-0.5 
-1.0 




z25 


none 
2.05^ 








i6q ; 

















-1.5 






r — \ 
j ^ 










\ -■ 


-2.0 








\ " 




i 




! " 
(1 : 


-2.5 
-J.O 






si 


\ " 




i 




11 J 
V 

'i : 


-0,5 
















-1.0 
















-1.5 




ii 

Ifi 
} 












-2.0 










- /' 






-2.5 
-3,0 






W 

M 
\l 

1..^ 




r J 

..i.i.. 








-1 


1 2 
V, fl[)^ km 5'^ 


3 4 




-1 


Q 


12 3 + 
[10^ km F"^] 



Fig. 15. — The same as Figure [M] but for model z25D. As with 
Model zl5D, a greater perturbation results in slightly enhances 
negative velocities, especially for ^^Ni. 
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Fig. 18. — Isotopic distribution as a function of mass for Model 
zl5D with different initial perturbations. A larger perturbation 
leads to more mixing. 
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Fig. 16. — Mass fraction of chemical isotopes as a function of 
radial velocity for model sl5A. Difference s be tween the models are 
larger than for Model s25A (see Figure [T7| but are not system- 
atic with perturbation, and become signiticant only below mass 
fractions of around 5% of the mass of ■'^"O and ^^Si. 
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Fig. 17. — Mass fraction of chemical isotopes as a function of 
radial velocity for Model s25A with different perturbations. There 
are only small differences, and the differences are not systematic 
with the amount of perturbation imposed on the initial model. 




Fig. 19. — Isotopic distribution as a function of mass for Model 
z25D with different initial perturbations. The distribution changes 
with the level of perturbation-larger perturbations lead to more 
mixing in a systematic way, and the difference between different 
perturbations is visible at the level of 50% mass fraction of ^He 
and -"^^O , much larger than in the more non-linear solar metallicity 
stars. 

Larger perturbations lead to a larger size scale for the ini- 
tial Ray leigh- Taylor fingers, which allows them to grow 
more quickly before the reverse shock passes them and 
the pressure gradient is no longer opposite the density 
gradient. In these models, the Ray leigh- Taylor insta- 
bilities cannot grow for many e— folding times, and so 
the initial scale of the perturbations matters. This can 
be seen in Figures [18] and [Toj which show the final dis- 
tribution of isotopes as a function of enclosed mass for 
different amounts of perturbation. The final amount of 
mixing is set by the size scale of the initial perturbation. 
The amount of mixing determines the distribution of iso- 
topes with velocity, as well, as shown in Figures 14] and 
[15I ''^^O and ^He are the most affected. The distance 
in mass space over which these isotopes are mixed varied 
by around 1 Mq for the perturbed models. 

3.3. Velocity Distribution 

The velocity distribution of chemical isotopes is shown 
in Figure 
in our so 
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^^Ti was mixed to relatively high velocities 
ar metallicity stars. The solar stars show 
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Fig. 20. — Isotopic distribution as a function of mass for Model 
sl5A with different initial perturbations. The isotopic distribution 
is not systematic with the amount of perturbation. -"^^O shows the 
greatest change, but only at mass fractions below 10%. 
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Fig. 21. — Isotope distribution as a function of mass for Model 
s25A with different initial perturbations. Where the distributions 
of isotopes are different, they are not systematic with perturba- 
tion. Differences between the runs are also small, and only become 
obvious for -"^^O , with oxygen concentrations are less than 10%. 

and "^He at high velocities of 4 x 10^ km s~^. These 
isotopes are also mixed all the way to the core of the 
solar metallicity stars, and some fraction of them reaches 
negative velocities of less than —0.5 x 10^ km s~^. In 
Model sl5A, where the most mixing takes place, we see 
^^O and ^^C mixed out to 2 x 10^ km s~^ in velocity 
space. ^^Ni does not reach the high velocities observed 
in 1987A. Model sl5A shows a peak in the ^^Ni velocity 
at ~ km s~^, with a tail extending to 0.7 x 10^ km 
s~^. Model s25A shows slightly higher velocities, with 
the ^^Nidistribution peaking at around 0.4 x 10^ km s~^ 
and reaching out to 0.9 x 10^ km s~^. These higher 
velocities for ^^Ni are probably due to smaller amounts 
of fallback in this star. 

The zero metallicity stars show lower velocities over- 
all. At the time they exploded as supernovae, these stars 
were more compact than solar composition stars. The 
supernova shock ran into a higher fraction of the total 
mass of the star at an earlier time. In neither zero metal 
model were the heavier isotopes mixed out to higher ve- 
locities than 1 X 10^ km s~^, but in Model zl5D the ve- 
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Fig. 22. — Mass fraction of chemical isotopes as a function of 
velocity for Models sl5A and s25A t = 9.26 days and zl5 and 
z25 at t = 2 and 4 hours, respectively. Negative velocities indicate 
material falling toward the point mass at the center of the star. 
Virtually no iron-peak elements escape from the primordial com- 
position stars, and their velocity distributions are smaller than the 
solar models. The peaks in the velocity distribution for ^^Ni for 
the solar stars are not high enough to match what was observed in 
1987A, but are consistent with previous studies of mixing in RSGs. 

locities of ^^Ni and ^^Ti had slight positive components, 
whereas in Model z25D the velocity distribution of these 
isotopes was almost entirely negative. No ^^Ni or ^^Ti 
e scape from Mode l z25D. 



[Kifonidis et al. 
wnich had a metal 



(2006) 



in their 1987A-type models, 
icity and mass intermediate to the 
stars studied in this paper, saw a well-mixed heavy el- 
ement core. The velocity profiles for isotopes from ^^O 
to the iron group were very similar. Because we did not 
see the iron core of our solar metallicity stars mixed to 
the degree seen in Kifonidis et al. (2006), our ^^Ni had a 
slower velocity distribution, while lighter isotopes, from 
^^Ti and ^^Si and lighter, were skewed towards higher 
velocities. 

The ^^Ni velocities observed in the zero-metallicity 
models were lower than those produced in earlier at- 
tempts at modeling SN 1987A, as well. The maximum 
iron -peak velocities obtained was on the order of 1300 km 



1991 



( Arnett et a l.|l989||Hachisu et al.ll990![ FYyxen et al. 
[Mueller e t al. 1991; Herant & Benz 199]]). In par- 
ticular, "Herant & Benz (1991) found that following the 
radioactive decay of ^°Ni and ^^Co increased the velocity 
of iron-peak elements slightly, but not enough to match 
observations of SN 1987A. Our iron-peak velocities were 
lower still. This was partly due to the large amount of 
fallback experienced by the primordial composition mod- 
els. Most of the ^^Nifell through the inner boundary of 
our zD-series simulations before it had time to power a 
nickel bubble like those seen in SPH simulations of SN 
198 7 A ( Herant fc Benz |1 991 ). Our progenitors were also 
more compact, with smaller helium cores, than the pro- 
genitor models for 1987A used in the above simulations. 

^^Ni was n ot mixed out to the the high velocities seen 
in SN1987A. |Utrobin| p004 ) report that ^^Ni was mixed 
X 10"^ km s~^. In no model did 
high velocities. Significant low- 



out to at least 2.5 
we see nickel at these 
order asymmetry in the explosion is probably necessary 
to reach ^^Ni velocities high enough to match observa- 
tions. 

In their 1994 paper on mixing in supernovae with red 
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supergiant progenitors, which employed progenitor mod- 
els similar to the sl5A and s25A progenitor models in this 
survey, Herant and Woosley found a similar distribution 
of isotopes in velocity space to the results presented here. 
They found a peak in the velocity distribution of ^^O at 
around 1.4 x lO^km s~^ for their 25 Mq star, we see one 
at around 1.1 x 10^ km s~^. We see oxygen mixed out 
to slight ly higher velocities i n our model si 5 A than are 



given m |.H.erant woosley| ( 19941). Their models show 
no negative velocities, which is probably a result of the 
SPH technique employed, which did not allow for accre- 
tion onto a sink particle at the center of the simulation. 

3.4. Yields 

The yields of our stars are greatly affected by the 
amount of fallback and degree of mixing induced by 
Ray leigh- Taylor instabilities. For models that experi- 
enced a great deal of fallback, as in case of the zero- 
metallicity models, mixing is of crucial importance in 
determining the final yields. Figure [3] shows the origi- 
nal distribution of isotopes in the four stars studied in 
this paper. The distribution of isotopes after all mixing 
has ceased is shown by Figure [4j As Figure |4] and Fig- 
ures |9] and [lO] show, layers interior to oxygen in the zero- 
metallicity stars have not mixed appreciably. A far larger 
portion of the solar metalicity stars has been mixed to- 
gether. The Ray leigh- Taylor instability has penetrated 
as far as the ^^Ni core of model sl5A, completely dis- 
rupting the exterior shells. This can be seen in Figure 
pT| Figure [12] also shows that Model s25A experienced 
more mixing than Model z25D (Figure 10), but not quite 
as much as Model sl5A. The lack of mixing into the in- 
terior layers of the zero-metallicity models means that 
very little of the silicon shell, and virtually none of the 
elements interior to this shell are mixed out. This, cou- 
pled with a large amount of infall, means that these Pop 
III stars expel almost no ^^Ni. 

To eliminate the possibility that the absence of 
Ray leigh- Taylor induced mixing at the edge of the ^^Ni 
shell in Models zl5D and z25D could have been due to 
mapping to two dimensions at too late a time we per- 
formed four simulations of the interior of Model z25D 
with an inner boundary of 1 x 10^ cm, an order of mag- 
nitude smaller than that used to simulate the star in the 
main simulations presented in this paper. These simula- 
tions were run at te n tim es the resolution of the simula- 
tions presented in § |2.4[ KEPLER models were mapped in 
10 and 40 seconds after bounce. To maintain consistency 
with the simulations that used a mapping from KEPLER 
at 100 s after bounce, we initially used a refiecting condi- 
tion at the inner boundary, which did not allow infall. At 
100 seconds after bounce for all of these high resolution 
models the inner boundary condition was changed to be 
zero- gradient. No Ray leigh- Taylor instabilities are seen 
to develop by 1x10^ s of simulation time, the point when 
the instabilities are clearly defined and growing along the 
C-O/He interface. Perhaps the perturbations from the 
Cartesian grid at this resolution are too small to effec- 
tively seed the instability. We perturbed the velocity on 
the grid with an / = 8 order Legendre polynomial with an 
amplitude 10% of the initial velocity. While this is proba- 
bly unphysical, it represents a limiting case: The absence 
of Ray leigh- Taylor induced mixing with a perturbation 
this extreme implies its absence in the larger, lower res- 



olution simulations is reasonable. The first 1 x 10^ s of 
the simulation were run at this high resolution. By this 
point the Ray leigh- Taylor instabilities are clearly visible 
at the C-O/He boundary, but there is no hint of the de- 
velopment of this instability at the interface betwe en the 
^^Ni core and the overlying layers. Kifonidis et al. (2006) 
speculated that the reason they saw mixing at the Si/O 
boundary when earlier studies did not was the "ad hoc" 
initialization of these explosions with a piston or thermal 
bomb, insufficient resolution, or differences in structure 
between their progenitor models and those employed in 
previous studies. Insufficient resolution does not appear 
to be the case here, and our extreme perturbations im- 
ply that seeding perturbations through neutrino-driven 
convection by starting the simulation earlier would have 
little effect. We can only conclude that this is a valid 
result arising from the structure of our progenitor star, 
or a result of the piston explosion mechanism. 

Because very little mixing occurs in the layers inte- 
rior to the mass cut of these zero metallicity stars, their 
yields were very sensitive to the amount of fallback that 
occurs. Figure [23] shows this sensitivity. When the mass 
of infalling mate rial was increa sed by 10% of the fiducial 
value taken from Zhang et al.| (|2008), almost none of the 
elements interior to the C-O shell escape. When the mass 
of infalling material was decreased by 10% of this fiducial 
value, more of the elements interior to oxygen escape and 
go on to enrich the surrounding gas. Figure [23] shows 
a comparison of the yields of our zero-metallicity stars 
with observations of the three most metal-poor stars yet 
found in the halo of the Milky Way. [X/H] values for 
all elements were obtained by diluting the yields of our 
model stars with enough pristine gas to fit observations 
of [C/H] and [0/H] to individual halo stars. The dotted 
and dashed lines show the yields resulting from increas- 
ing or decreasing the amount of fallback by 10% of its 
fiducial value, respectively. Model zl5D is more sensi- 
tive than Model z25D to decreasing the mass cut. The 
model yields appear much closer to fitting the values of 
elements heavier than oxygen observed in the metal-poor 
halo stars. 

The yields of our Pop III core collapse supernovae re- 
produce the extreme overabundance of [C,0/Fe] to the 
point where Fe is underproduced even when compared to 
the already low observations of [C,0/Fe] in HMP stars. 
Carbon is slightly underproduced and nitrogen is under- 
produced by 2-3 dex relative to oxygen, as is the iron 
peak. The iron peak elements are underproduced be- 
cause they do not mix sufficiently with the lighter el- 
ements in the star and cannot escape falling back onto 
the remnant at the center. A reduction in the infall mass 
brings zl5D closer to reproducing observations. Two di- 
mensional versions of the one dimensional simulations 



presented in th e high er energy exp l osions of |Heger fc 



Woosleyl ( |20Q8| and N omoto et ah] ( |2007D might eject 
more of the ^"Ni produ ced in the more ene rgetic explo- 
sions. Rotating models ( Hirschi et al.||2008 ) create more 
primary nitrogen, which can lead to an increase in the 
rate of CNO burning at the base of the hydrogen shell, 
causing some models die as larger red supergiants rather 
than a compact blue supergiants. In this case, Rayleigh- 
Taylor mixing would play out in a similar way to the solar 
models presented in this paper, and more iron would be 



ejected by these stars. Recent simulations (Scheck et al. 
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Fig. 23. — Model results matched to carbon and oxygen abun- 
dances found in individual HMP halo stars. Model zl5D, with 
less fallback, expels some iron peak elements. Nitrogen is always 
under-produced, which is not a result of insufficient mixing during 
the supernova explosion but rather the presupernova evolution of 
these models. The solid line is the yield calculated using the rem- 
nant masses from Zhang et al. (2008). The dotted and dashed lines 
correspond to the amount ot mfall found in Zhang et al. ( 2008 ) de- 
creased or increased by 10%, respectively. The sensitivity of the 
yields to fallback is especially apparent for Model zl5D, the model 
with the least mixing. 



2QQ4[ |2QQ6| [Burrows et 'aI]|2QQ7a|b|cl ) point out that the 
supernova explosion mechanism is probably inherently 
multidimensional and asymmetric. Asymmetry in the 
explosion, whether in the form of a jet or a perturbation 
described by Legendre polynomials of order of / = 1 or 
/ = 2, might also mix more of the nickel core out of the 
star, brin^in^ the mode l s clos er to reproducing observa- 



tions. Venn & Lambert (2008) have suggested that HMP 
stars may be "chemically peculiar" stars, in which low 
iron abundance is caused by separation of gas and dust 
beyond the stellar surface, followed by accretion of dust- 
depleted gas. If this is the case-and the authors note 
that a definitive answer requires additional information- 
the stars' true metallicity is closer to [X/H] ^ —2 rather 
than -5. 

3.5. Visibility 

The supernova light curve is affected by the amount 
of ^^Ni in the center of the star that falls back onto the 
black hole at the center of the explosion. The models 
for the first supernovae presented in this work are in- 
trinsically dimmer than corresponding supernovae aris- 
ing from stars of solar composition provided they ex- 
plode with the same amount of energy. In our models 
of primordial composition supernovae, all or nearly all of 
the ^^Ni synthesized in the supernova falls back onto the 
remnant left behind at the center of the explosion. En- 
ergy from the radioactive decay of ^^Ni powers the tail 
of core-collapse supernova light curves. When the en- 
ergy released in its radioactive decay to ^^Fe is no longer 
observable, the supernova light curves will loose their ra- 
dioactive tails, making them briefer and dimmer than 



ordinary core-collapse supernova light curves. 

4. CONCLUSIONS 

The presupernova structure of a star is determined 
largely by its initial mass and by the initial composi- 
tion of the gas from which it formed. The symme- 
try and energy of the explosion, along with the presu- 
pernova structure, infiuence where and to what extent 
Ray leigh- Taylor instabilities will grow, as well as how 
much mass will fall back onto the remnant at the center. 
The non-rotating zero metallicity models studied are far 
more compact than solar-composition models of the same 
mass, in part because CNO burning proceeds at higher 
temperatures and densities. CNO burning is responsible 
for energy production during the main sequence for all 
stars at the masses studied here, but in metal poor stars 
CNO burning proceeds at higher temperatures and den- 
sities. For zero-metalicity stars, the star must first con- 
tract to a temperature of 10^ K, hot enough to initiate 
helium burning. This helium burning produces a small 
amount of carbon, which is enough to act as a catalyst to 
enable hot CNO burning to proceed. In addition, non- 
rotating stars with a metallicity Z below 10~^ will never 
reach the red giant branch, since they end helium burning 
with effective temperatures above 10^. Below this tem- 
perature, the opacity is large enough that the star will 
expand toward the red giant branch. The more compact 
structure of these stars causes their reverse shocks to 
propagate more quickly to the origin than those in solar 
stars. Larger remnants are left behind in the more com- 
pact stars because the rate at which mass acc retes onto 
the stellar remnant is higher, as predicted by [Chevalier 
(il989t) and shown in the ID simulations of Zhang et al. 



The time scale over which the Ray leigh- Taylor insta- 
bilities can develop is also set by the reverse shock. For 
the case of the compact primordial composition progeni- 
tors modeled here, the Ray leigh- Taylor instabilities have 
little time to develop. This means that a smaller por- 
tion of the isotopic layers of the star will be mixed. The 
Ray leigh- Taylor instabilities do not have time to become 
fully nonlinear in our simulations, so the scale of the in- 
stability as well as the degree of mixing is set by the scale 
of the initial seed perturbations. In the case of the solar- 
composition progenitor models, the Ray leigh- Taylor in- 
stability became fully nonlinear and the size and shape 
of the initial perturbation was no longer apparent at late 
times. A smaller region of the primordial-composition 
stars is unstable, compared to solar-composition stars, 
which also contributes to the reduced mixing we see in 
our zero- metal models. The small amount of mixing 
experienced by zero-metallicity stars means that their 
yields are very sensitive to the amount of fallback, which 
may depend somewhat on the time the models were 
mapped from KEPLER to FLASH. 
Unlike the recent simulations of Ray lei^h- Taylo r mix- 



ing in compact blue supergiants by Kifonidis et al. (2003 
2006), we do not see mixing at the silicon-oxygen shel 



interface. The absence of instability at that interface, as 
well as at the nickel-silicon shell interface in zD stars is 
robust. Mixing was not observed even after mapping the 
models to FLASH at earlier times, running the models at 
significantly higher resolution, and adding large pertur- 
bations of 10% to the velocity. Differences between the 
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progenitor models used in Kifonidis et al. (2003, 2006) 
and our progenitor models, or structural differences aris- 
ing from their early modeling and different explosion 
mechanism may be responsible. The degree of mixing 
we observe in our models is similar to that of earlier 
studies of compact blue supergiant progenitor models for 
SN 1987A. The degree of mixing and final velocities of 
isotopes we observe in our s olar composition m odels is 
comparable to that found in lHerant fc Woosleyj ([l994). 
Material from the silicon and nickel layers is mixed out 
into the outer layers of the solar metallicity stars, and 
these isotopes are mixed to velocities of ~ 1 — 2 x 10^ 



km s vel ocities comparable to those found in Herant 
Woosley| ([T994). 



The zero-metallicity models show a dramatic overpro- 
duction in carbon and oxygen relative to iron. This is 
the trend observed in the most metal-poor stars in the 
universe, but the degree of mixing in the present work is 
so small that not enough iron escapes to reproduce the 
already large [C+O/Fe] values observed in these stars. 
Decreasing the amount of fallback onto the remnant at 
the center through a more energetic explosion might suc- 
ceed in driving more iron-peak isotopes out of the star, 
but would not result in a higher amount of nitrogen. Our 
models also cannot reproduce the [C/N] or [0/N] ratios 
observed in HMP stars. Zero-metallicity stars with sig- 
nificant rotation might produce both more primary ni- 
trogen as well as more iron. Rotating models experience 
more shear mixing, which can drag primary carbon and 
oxygen from the ^He burning convective core to the 
burning shell, where it becomes nitrogen, increasing the 
amount of burning in the shell and making that shell 
convective. The zero-metallicity st ars may then die a s 
a red, rather than blue, superg iant dHirschi et a'L]|2008D . 



If this is the case, then rotating zero-metallicity models 
might look very much like the solar-metallicity models 
presented in this paper, which would more closely repro- 
duce abundances observed in HMP stars. A higher SN 
explosion energy might be another way to expel more 
iron from these stars. 

Future work will include modeling primordial composi- 
tion models wit h some degree of rotation, as recent work 
( Hirschi et al.|2 008) indicates that they should have more 
primary nitrogen, as well as less compact structure and 
a larger ^He core. Modeling rotating stars may result in 
drastically different yields for our zero metallicity stars, 
allowing more iron to escape and increasing the [C/N] 
ratio closer to what is observed in HMP stars. We will 
explore higher explosion energies and asymmetric explo- 
sions and how they affect the mount of mixing. 
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